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ABSTRACT 

Context. The existence of large-scale and long-lived 2D vortices in accretion discs has been debated for more than 
a decade. They appear spontaneously in several 2D disc simulations and they are known to accelerate planetesimal 
formation through a dust trapping process. In some cases, these vortices may even lead to an efficient way to transport 
angular momentum in protoplanetary discs when MHD instabilities are inoperative. However, the issue of the stability 
of these structures to the imposition of 3D disturbances is still not fully understood, and it casts doubts on their long 
term survival 

Aims. We present new results on the 3D stability of elliptical vortices embedded in accretion discs, based on a linear 
analysis and several non-linear simulations. 

Methods. We introduce a simple steady 2D vortex model which is a non-linear solution of the equations of motion, 
and we show that its core is made of elliptical streamlines. We then derive the linearised equations governing the 3D 
perturbations in the core of this vortex, and we show that they can be reduced to a Floquet problem. We solve this 
problem numerically in the astrophysical regime, including a simplified model to take into account vertical stratification 
effects. We present several analytical limits for which the mechanism responsible for instability can be explained. Finally, 
we compare the results of the linear analysis to some high resolution numerical simulations obtained with spectral and 
finite difference methods. A discussion is provided, emphasising the astrophysical consequences of our findings for the 
dynamics of vortices. 

Results. We show that most anticyclonic vortices are unstable due to a resonance between the turnover time and the 
local epicyclic oscillation period. A small linearly stable domain is found for vortex cores with an aspect-ratio of around 
5. However, our simulations show that it is only the vortex core that is stable, with the instability still appearing on 
the vortex boundary. In addition, we find numerically that results obtained under the assumption of incompressibility 
are not affected by the introduction of a moderate compressibility. Finally, we show that a strong vertical stratification 
does not create any additional stable domain of aspect ratio, but it significantly reduces growth rates for relatively 
weak (and therefore elongated) vortices. 

Conclusions. Elliptical vortices are always unstable, whatever the horizontal or vertical aspect-ratio is. The instability 
can however be weak and is often found at small scales, making it difficult to detect in low-order finite-difference 
simulations. 

Key words, accretion, accretion disks - instabilities - hydrodynamics 
1. Introduction rotational instability (jBalbus fc Hawlcv 1998) doesn't op 



■ erate, such as in dead zones (|Gammielll996| ). Second, they 



The existence of 2D lo ng-lived vortices in ac cretion discs are a vei T efficient way to accelerate t he planetesimal 
was first proposed by Ivon Weizsackerl (| 19441) in an out- form ation process in protoplanetary discs flJohansen et al.| 
moded model of planet formatio n. This idea was re- 12Q0J). However, the stability of these vortices when small 
vived bv lBarge fc Sommerial (|l995l ) to accelerate planetes- 3D disturbances are imposed is largely unknown, 
imals formation by a dust trapping process. This kind In the astrophysics comm unity, this issue ha s been in- 
of vortex is often observed in 2D simulations of discs vestigated mainly numerically. IShen et all (|2006l ) examined 
(see e.g. iGodon fc Livid Il999t lUmurhan fc Reeevl l2004t the formation of 2D vortices starting from 2D turbulence in 
I Johnson fc Gammidl2005bl : lBodo et al.ll2007l ). since 2D tur- fully compressible simulations. According to their results, a 
bulence is known to generate a n inverse casca de of en- small 3D noise added to their initialy 2D configuration de- 
ergy leading to large 2D vortices (|Onsagerl [l949). Vortices stroys the coherent v ortices in a few orbits , relaxi ng the flow 
may also be gene rated by 2D instabili ties such as Rossby to its laminar state. iBarranco fc Marcus! (|2005f) also corn- 
wave in stabilities ( Lovelace et al. 19 991) or baroclinic insta - puted the evolution of 3D vortices usin g an anelastic cod e 
bilities (Kla hr~fc Bodenheimerll2~003l : iPetersen et al.l [2007). incorporating vertical stratification. As IShen et alJ ((2006), 
although the latter is stil l a m atter of ongoing debate they found that midplane vortices were destroyed by 3D 
(1 Johnson fc Gammiel[2005all2006f) . These vortices may play perturbations. However, they also showed that off-midplane 
at least two important roles regarding accretion disc dy- vortices could survive for several hundreds of orbits, leading 
namics. First, they could lead to an efficient angular mo- to the possibility of a stabilizing effect due to the stratifi- 
mentum transport process in regions in which the magneto- cation 
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It is often assumed that these vortices are unstable 
because of the elliptical instability. The elliptical insta- 
bility is a parametric instability appearing when a mul- 
tiple of the vortex turnover frequency matches an iner- 
tial wave frequency, leading to a positive resonance. It is 
observed when the backgound flow follows closed stream- 
lines, and being localized on individual streamlines is a 
local instability (in particular it doesn't need to involve 
the vortex bo undaries). This instab ility was first found 
numerically bv iPierrehumbertl (|l986l ) an d desc r ibed u sing 
ICraik fc Criminald (119861 ) solutions by iBavlvl (1 1986ft for 
pure e lliptical flows. The rotating case was studied bv lCraikl 
(|l989ft . who showed that anticyclonic elliptical flows can be 
stabl e for some rota tion rates. Interested readers may con- 
sult E erswc 3 (120021) for a more extensive discussion of the 
elliptical instability and its development in fluid mechanics. 

In the present paper, we investigate the elliptical in- 
stability in the context of accretion disc vortices. We first 
present a steady 2D vortex model, which is a non-linear 
solution of the local disc equations. We then present the 
linearised equations governing 3D perturbations inside the 
vortex. A criterion for the instability is derived from these 
equations and a physical understanding of the mechanism 
responsible for the instability is provided. We briefly ex- 
tend these results to a simplified stratified case, and we 
compare our findings to fully non-linear simulations of ac- 
cretion discs vortices. Finally, we provide a discussion and 
a comparison with previous work. 



2. Local model of an elliptical vortex. 

2.1. Shearing-sheet model 

In the following, we will assume a local model for the accre- 
tion disc, following the shearing-sheet a pproximation. The 
reader may consultlHawlev et al.l (|l995ft . lBalbusl |2003) and 
iRegev fc UmurhanP 20081 ) for an extensive discussion of the 
properties and limitations of this model. As a simplification, 
we will assume the flow is incom pressible, consistently wit h 
the small shearing box model (|Regev fc Umurhanl [2008) . 
The shearing box equations are found by considering a 
Cartesian box centred at r = Ro, rotating with the disc 
at angular velocity £1 = £l(i?o). We define Ro<fi — + x and 
r — i?o - * — y for consistency with the standard n otation 
for plane Couette flows (e.g. IDrazin fc Reidl [l981h . Note 
that this definition differs from the stan dard notation used 
in shearing boxes (jHawlev et al.l [l995h with x — > —ysB, 
y — * £sb and z — ► zge- In this rotating frame, one obtains 
the following set of governing equations 



d t u + V • (u (£) u) 
V • u 



-VII - 20 x u + 2flSye y , 



0. 



(1) 

(2) 



In these equations, we have defined the mean shear S = 
—rd r Cl, which is set to S = (3/2)17 for a Keplerian disc. 
The generalised pressure II = P/po is calculated solving a 
Poisson equation with the incompressibility condition. One 
can check easely that the velocity field u — Sye x is a steady 
solution of these equations. 



subsection. Since we are looking for 2D solutions in the 
(x, y) plane, we can omit the Coriolis force as it will only 
change the pressure distribution. We define this vortex by 
an elliptical patch of constant vorticity uj t = —S + u v (the 
"core") where — S is the background flow vorticity and ui v 
is the vorticity of the vortex itself. Outside of this core, the 
vorticity is a ssumed to b e u>t — ~S, extending to infinity. 
According to lKidal (|l98lt ) (Eq. 2.9), such a vortex is steady 
if the semi-major axis is aligned with x and if its vorticity 
satisfies 



~S 



x + i 



-If 

x v x - 1 



(3) 



where we have defined the vortex aspect-ratio x — a l°i a 
and b being respectively the vortex semi-major and semi- 
minor axis. One deduces from this result that only anticy- 
clonic vortices (CI and lo v having opposite sign) are stable 
in Keplerian flows. Since this solution is steady, no stream- 
line goes through the core boundaries, and the streamlines 
inside the core have to be elliptical, with the same aspect- 
ratio as the vortex core. 

Thanks to this property, we can write the velocity field 
in the vortex core, assuming it's centered on x — y = 0, as 



l 



-xv, 



Uy (x-i)x 



1 

— x. 



This solution can be written in the simpler form 



— S A.ij xj , 



defining 



(4) 
(5) 

(6) 
(7) 



2.3. Explicit solution 

A complete solution for the velocity field can be found 
defining the streamfunction ip(x, y) so that u x = —d y ip and 




d x ip. This streamfunction satisfies: 



— S + lo u inside the core, 
—S outside. 



(8) 



One solves these equations in elliptical coordinates, defining 

(p, v) by 



x = f cosh(/i) cos(^), 
V = /sinh(/i)sin(^), 



(9) 
(10) 



with / = a^(x 2 - 1)/X 2 
boundary is found at p 



In these coordinates, the core 
= po with tanh(/io) = x^ 1 ■ 
Requiring that ip and d^ip are continuous at po, one finds 
the following expressions for ip: 



2.2. The lKidl Hl98l) solution 



We want to study the stability of a steady elliptical vortex 
embedded in the sheared flow described in the previous 



ipi 



Sf 



2(X- 



— (x 1 cosh 2 (p) cos 2 (v) 



+X sinh 2 (/i) sin 2 (r/ 



(11) 
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Fig. 1. Vortex solution streamlines in a sheared flow with 
X = 3. As expected, inside the vortex (delimited by the 
blue bold line) the streamlines are elliptical, with the same 
aspect ratio as the vortex itself. 



i> 



sf 



4( X -1) 2 



1 + 2(/i - /i ) 
+2(x- l) 2 sinh 2 (^) sin 2 (j/) 



x-i 



X ■ 



-exp[-2(^-^ )]cos(2i/) , (12 



where the subscripts i and o stand for inside and outside 
the core. 

As expected, the inner solution reproduces the ellipti- 
cal flow described in the previous subsection (Eqns. @]— [5]). 
In the outer solution, one recognises the background shear 
(3rd term), which dominates for large fi. Moreover, this 
solution exhibits a linear term in /x corresponding to the 
long-range perturbation of the vortex. In cartesian coordi- 
nates, this linear dependance translates into a logarithmic 
tail for ip(x,y), showing that the vortex presence can be 
felt far from the core. As we will see in section [5] below, this 
property leads to numerical artefacts when one tries to fit 
this solution in a finite-size box. For completeness, we show 
in Fig. [1] the streamlines obtained from solution (TTTj) — (|12p . 
As expected, the streamlines inside the core are ellipses of 
constant aspect ratio x- Outside the core, one still finds 
closed streamlines but with a much more elongated struc- 
ture. 



2.4. Perturbation equations 

In the following, we concentrate on the evolution of pertur- 
bations in the vortex core, assuming the latter is infinite. 
This corresponds to a limit in which the perturbations are 
small compared to the typical horizontal size of the core. 
To study the evolution of 3D perturbations in a 2D ellip- 
tical flow, we write the total velocity field as u = u° + v 
where v i s suppo s ed to have an infinitely small amplit ude. 
Following EiM^ (|1880T) and lCraik fc Criminald (|1986D . we 
use a time explicit Fourier decomposition for the perturba- 
tion v = v(t) exp(ik(t) ■ x). Plugging this solution in (fl]) 
leads at first order to: 



it + ix k Vi(k k + SkjAjk) = /"/.-, 1 1 Sr,.\,, 

^^ij k ^k 

kjVj = 



(13) 
(14) 



where we have included the tidal term of fTJ) in II. To satisfy 
(fHi|) for any x k , one has to solve: 



k k + SkjAjk = 0, 



(15) 



which leads to: 

fe(i) = fco^sin(0) cos [^(t)]e a 

-Xsin(0) sin [<f>(t)]e y 

+ cos(8)e z 



(16) 



where ko, and to are integration constants and <p(t) = 
S/(x ~ — *o) is the turnover angle. Therefore, the per- 
turbations will have a rotating wavevector with a turnover 
time equal to the vortex turnover time T = 2ir(x — !)/£• 
Because of the incompressibility condition, 9 = will corre- 
spond to horizontal motion perturbations whereas 9 = tt/2 
will imply vertical ones. We also note that the rotating 
wavevector will involve smaller structures in the y direc- 
tion. According to this result, the horizontal aspect ratio of 
the perturbation wavevector is equal to that of the vortex 

(x). 

We can then simplify (|13[) eliminating the pressure, 
which leads to the final system 



(^ + ^ -%) 



hi kj 



(17) 



where A = (x~l)A and Rj m = (%- l)e^ m O;/5. Plugging 
solution (flB) in this equation leads t o a Fl oquet problem 
for v, as already pointed out bv lBavlvi (|1986f ). Interestingly, 
the Floquet problem doesn't depend on the norm of the 
wavevector k = ko, but just on its direction. Therefore, this 
problem is scale- independant, at least in the inviscid limit. 
The solution to this problem is known to be a superposition 
of Floquet modes written 



v = cxp( 7 t)/[<Xi)], 



(18) 



where / is periodic with period 2n. To determine the 
Floquet exponents 7, we compute the eigenvalues exp( 7 T) 
of the matrix M(2tt) where M(<fi) satisfies the generalised 
equation 



dM,, 



(CW~~ %)^m + 2(^-^)^ jro ] M mn , (19) 



with the initial condition 

M ij (0)=5 ij . (20) 

A numerical approach is required to solve this problem in 
most cases. However, some limits can be understood using 
analytical approaches, as we will see in the next section. 



3. The elliptical instability 

3.1. Horizontal instability 

It is possible to derive an analytical criterion for the insta- 
bility in the limit k = k z e z . In this limit, v z — and (flT]) 
is reduced to 



dv x 
~d4 
dv„ 



X 



20, 

T (x-D 



1 20 



(21) 
(22) 



This system describes horizontal epicyclic oscillations with 
frequency 



k 2 = S 2 { R 



X 



R 



1 



x(x- 1) 



(23) 
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having defined the rotation number as R = 20/ S. In the 
limit x ~ > 00 where the vortex is infinitely weak, we find 
the classical epicyclic frequency k 2 — 20(20 — S), which 
is equal to the Keplerian frequency in discs. This epicyclic 
frequency is imaginary when 



which direcrly leads to 
that if 



To show this we first note 



and /3=(^ + 20 



(30) 



y/T+4/R 



<X< 



R 



R-l 



(24) 



If a columnar vortex is in this regime, its horizontal layers 
will tend to drift exponentially in the (x, y) plane, leading to 
the destruction of the vertical coherence of the structure. In 
accretion disc vortices, this regime is found for 3/2 < x < 4. 

3.2. A generalization 

We note that the above discussion may be generalized to ap- 
ply to a wider class of steady flows (u x , u y ). The linearized 
equations of motion are 



where II' is the pressure perturbation. The operator 



^ d 



(26) 



denotes the convective derivative. We now assume pertur- 
bations are local in z so that we may assume that the z 
dependence is through a factor exp(ik z z), where k z is very 
large in the sense that the length scale fc~ x is smaller than 
any other in the problem. From the incompressibility con- 
dition V • v = 0, and the z component of ((25)) we conclude 
that II' = 0{k~ 2 ). Thus it may be dropped from the x and 
y components of (|25[) which then give the pair 



Vv x 



du° x 
dx 



du° y 

'~dy- 



'dul 



dy 



dx 



^ -20 



20 



(27) 



(28) 



Replacing T> by 4?, with the understanding that the time 
derivative is to be taken on a fixed streamline, we see that 
in general the system (|27)) - (l28l) leads to a second order or- 
dinary differential equation with periodic coefficients, the 
period being the time to circulate round the chosen stream- 
line. But note that for the problem on hand the coefficients 
are constant because the unperturbed velocity is a linear 
function of the coordinates. In more general cases one has 
a Floquet problem of the type described above for every 
streamline indicating that unstable modes are indeed local- 
ized on streamlines. In order to connect with (|24[) we note 
that one can prove that if everywhere on a chosen stream- 
line 



dy 



20 



du° y 
dx 



20 > 0, 



(29) 



there will be an instability. In fact for the problem on hand 
this is precisely equivalent to the above condition that k 2 < 



then both a and (3 must be either positive or negative. 
Without loss of generality we assume both to be positive 
(the case when they are both negative may be recovered by 
setting v x — > —v x ). Then from equations ([27|) -([28 |l we may 
obtain 



d(v x v y ) _ 
Jt 

or setting q y 
d(v x q v ) 



av y - (3v x 



dt 



aq y + (3v x 



(31) 



(32) 



Thus if 

Qmin and f3 m i n are the minimum values of ot and [3 
respectively, we have 



( 25 ) d(v x q y ) 



< 



dt 

or equivalently 



(33) 



d(v q ) ^ 

Thus it follows that v x q y grows faster than exponentially 
with growth rate 2^J 'a m i n (3 m i n which means that there is a 
Floquet exponent corresponding to exponential growth. 



3.3. Numerical results 

In the general case (9 > and x > 1), the Floquet prob- 
lem (|T71) can't be solved analytically. We therefore use the 
following numerical approach. We first evolve the M ma- 
trix following Eq. [19] for a given \ an d with a 4th order 
Runge-Kutta-Fehlberg algorithm, keeping the error down 
to 10 -10 . The eigenvalues are then computed from M(2tt) 
using the GNU scientific library routines. This procedure 
is repeated for 1000 values of x and 1000 values of 9 to get 
a 2D representation of the regions for which the instability 
is found. 

To test our nume r ical re sults, we have first reproduced 
the results of iBavlvl §986) for a vortex in a non rotat- 
ing sheared flow. One find in this case that vortices are 
unconditionally unstable to fully 3D disturbances (Fig. [3]). 
Moreover, in the limit ^ — > 1, only modes with 9 = 60° are 
found to be unstable, consistently with Bayly's findings. 
Finally, we note that the growth rate for the instability 
weakens as we elongate the vortex. This finding, although 
apparently opposite to the classical result for the ellipti- 
cal instability, is due to the fact that our growth rate is 
measured in inverse shear time and not in inverse vortex 
turnover time. 

We have used our numerical method to solve the 
Floquet problem in the Keplerian case (Fig. [3]). In this 
case, we observe essentially a two bands structure. The first 
band (low-x band) is located between x = 1 an d x = 4. 
It is associated with a short growth time (typically of the 
order of one shear time) , and it includes the horizontal in- 
stability described in section 13.11 for 9 = 0. In the limit 
X —> 1 (infinitely strong vortex), we find the same result 
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5 10 15 

X 

Fig. 2. Growth rate (in shear time) of the 3D instability 
deduced numerically from the Floquet problem (fT7|) in the 
non rotating case. This result i s ident i cal to the original 
elliptical instability described bv lBavlvl |l986) and allow us 
to test our numerical approach. The colour scale represent 
the logarithm of the growth rate. 



as in the non rotating case (instability for 9 = 60°), which 
was to be expected since the vortex turnover time is much 
smaller than the rotation time in this limit. The second 
band (high-x band) is found for x > 6 and is much weaker 
since it involves growth rate of the order of 10 _2 5. This 
band gets wider and weaker as we go to large aspect ratios 
(and weak vortices), with a growth rate maximum found for 
X — 11. We have tried to get a larger resolution in the re- 
gion x ~ 5.5 — 6 and it seems that this band reaches 9 = 0° 
for x ~ 5.9, but the growth rates involved are very small. 
Note that another narrower band is observed in the regime 
of high x an d smaller 6 in Fig. [3] As we will see in the 
following, this band is due to higher order resonance and is 
therefore weaker than the bands described previously. 

To get a simpler representation of the elliptical insta- 
bility in the Keplerian case, we have plotted the maximum 
growth rate as a function of x hi Fig. [U We observe the 
same basic features as on the colormaps with 2 bands. We 
also observe that the high-x band decays as x g ets larger. 
This property is to be expected since the limit x ~ > 00 
corresponds to a pure shear flow, which is known to be lin- 
early stable. However, for a finite (but large) x, vortices 
are always uns table, contrary to claims in the literature 
(lLithwicldl200l . 

Therefore, it seems that no elliptical instability is found 
for 4 < x < 5.9. However, we recall that up to now we 
have only considered perturbations interior to the vortex 
core. But as we will see later, instabilities may also appear 
outside the vortex core. 

3.4. Physical interpretation 

It may be noted that the system ([XT]) can be writte n as 
a single Hill's equation (see for example IWaleffd ["l990l) ■ In 
this case, one finds that the periodic excitation frequency 




2 4 6 10 20 40 80 



Fig. 3. Growth rate of the 3D instability deduced numeri- 
cally from the Floquet problem (TT7|) in the Keplerian case 
for large x- One observes the vertical modes for 3/2 < x < 4 
and full 3D modes for x < 3/2 and x > 6 up to at least 
X = 100. A weak secondary band is found under the pri- 
mary band for large x- 




10 1 
X 

Fig. 4. Maximum growth rate in the Keplerian case (in 
shear time). We remark the low-x band with large growth 
rates (x < 4), the high-x band with smaller growth rates 
(x > 6) and a stable region in between. The growth rate 
decays for large x when the vortex gets weaker, as expected. 



is twice the vortex frequency S/(x — 1), since the time- 
dependent wave numbers k x and k y always appear quadrat - 
ically. In the limit 9 — > for which the excitation is small, 
one would expect a series of instability bands. By analogy 
with Mathieu's equation, one should find such a band when 
the local epicyclic frequency n is an harmonic of the vortex 
frequency: 

k = n -, (35) 

x- 1 

with n e JV. In the Keplerian case, one would therefore 
expect instability bands arising from 9 = for xi — 4.65, 
Xi — 5.89, X3 — 7.28. . . However, according to Fig. [3l the 
first unstable band x = 4.65 (n = 1) doesn't appear. This 
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peculiar property can be understood writing the epicyclic 
frequency as 



K 2 



RT 



1 



(x-i) 2 



(36) 



where r = u>t/S + R is the absolute vorticity of the back- 
ground vortex. According to this expression, the n = 1 
resonance condition is equivalent to the constraint T = 0. 
Interestingly the absolute vorticity also appears explicitly 
in the linearized vertical vorticity equation: 



d t uj z + Mo • Vw 2 = STd z v z 



(37) 



where lu z = d x v y 



d v v x is the vertical vorticity of the per- 
turbation. According to this equation, u z is a conserved 
quantity if the background absolute vorticity T is zero. We 
conclude from this argument that the n — 1 resonance con- 
dition cannot lead to an instability, since oj z could not be 
conserved in that case. Therefore, the first unstable band 
appears for n = 2 (x — 5.89), which is consistent with 
our numerical result. As stated before, higher order bands 
(n > 2) are also observed in the computation (Fig. [3]), but 
they are much weaker, and one can't check the origin of 
these bands on the 8 = axis with enough precision to 
compare with eq. ([33)1 . 



4. Stratified case 

4.1. Model and Equations 



iBarranco fc Marcus! (|2005f ) showed using anelastic numeri- 
cal simulations that off-midplane vortices were able to sur- 
vive for several hundred orbits whereas midplane vortices 
were destroyed rapidly. This suggests that stratification 
might be a way to stabilise vortices and suppress the el- 
liptical instability. To study this problem in a simple con- 
figuration, we have considered the case of a 3D vortex cen- 
tered at z = zo, i.e. above the midplane. We then follow the 
evolution of perturbations with k 3> Zq 1 inside the vortex. 
In first approximation, the evolution of these perturbations 
can be described using the Boussinesq approximation which 
can be written 



d t u + V • (u i 



V • «) 
V • u 



vn- 

-NCe, 
-Nu z , 



2ft x u + 2ClSye v 



0. 



(38) 

(39) 
(40) 



where £ is the potential temperature and ./V is the Brunt - 
Valsala frequency. In the following, we will assume a stable 
vertical stratification, i.e. a real Brunt- Vaisala frequency. 
The vortex solution found previously is still valid in the 
Boussinesq approximation. Therefore, we can follow the 
same procedure as the one used in the non stratified case 
to get 



dvi 



d( 



( 2 k^ k j 
N( X - 



Sij j Aj m 

-1) 



(i^-%)c<^ 



N( X -l) 



(41) 



(42) 



These equations associated with the solution (|16|) lead to 
a stratified Floquet problem. As previously, this problem 




X 

Fig. 5. Growth rate of the 3D instability deduced numeri- 
cally from the Floquet problem (|4"Tj) in the Keplerian case 
with N/S = 1.0. One observe once again the vertical modes 
for 3/2 < x < 4 and many weak bands involving fully 3D 
perturbations. The colour scale represent the logarithm of 
the growth rate. 



can be solved by diagonalizing a 4 x 4 evolution matrix M, 
following the same procedure as in the non stratified case. 
The numerical method used to solve these equations is the 
same as in the previous section, except we have added an 
extra dimension in the solver to handle the evolution of the 
potential temperature. 

4.2. Results 

We plot the growth rate of the instability as a function of 
X and 8 for a moderately stratified vortex N = S in Fig. [5l 
We note that the influence of the stratification is strong in 
the domain in which the vortex is weak (% > 5), and signif- 
icant differences are observed when comparing with Fig. [3J 
First, we remark that the stable region observed in the non 
stratified case (4 < x < 5-9) is not present in the N = S 
stratified case. This result is due to the presence of high 
frequency buoyancy modes which can match the turnover 
frequency when the inertial modes cannot. Moreover, the 
high-% band has been replaced by a series of weaker bands. 
These bands can be understood as reminiscences of the res- 
onance condition (f35|) . since in the limit 8 — > the effect of 
stratification disappears. 

When we compute the maximum growth rate as a func- 
tion of x f° r various stratification frequencies (Fig. [6]), we 
note that the growth rate at large x decreases significantly 
with stronger stratification. In this sense, the stratification 
tends to weaken the elliptical instability, but does not sup- 
press it. 

5. Non-linear simulations 

5.1. Resolving the elliptical instability 

As pointed out in the previous section, the elliptical in- 
stability is always present for highly elongated vortices 
(x > 6), although it is quite weak. We have also found that 
no elliptical instability was observed for 4 < x < 5.9 when 
stratification was omitted. To check these results, one has 
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Fig. 6. Maximum growth rate in the Keplerian case (in 
shear time) for various stratification frequencies. As the 
stratification gets stronger, the growth rate decreases for 
weak vortices (x > 6) . We also note that for a strong enough 
stratification, the region without instability (4 < \ < 5.9) 
disappears. Finally, the horizontal instability described in 
the previous section (3/2 < \ < 4, 9 — 0) is not affected 
by the stratification, as expected. 



to perform non-linear simulations in which this instability 
is resolved. A condition to get this instability in a numeri- 
cal simulation may be found assuming we have a 2D vortex 
embedded in a disc with a height H , the vortex core being 
larger than H . In the limit \ ^ 6, we find the elliptical 
instability for 9 ~ 7r/3 (Fig. [3]). According to (I16|) . unstable 
modes will have in this case fc™ ax ~ k z and fc™ ax ~ X^z- 
Moreover, k z > 2ir/H since 3D perturbations have to fit 
vertically inside the disc. If wc assume wc need n points 
to resolve one wavelength (n > 2), then the longest wave- 
length modes unstable to the elliptical instability will be 
resolved if the x resolution is n points per scale height and 
the y resolution is n\ points per scale height. This last con- 
dition is not often satisfied in 3D global simulations, and 
may explain why long-lived and stable vortices are some- 
times observed. Note that the value n may be quite high 
when using low order finite difference methods (one might 
expect n ~ 10) since one wants the (numerical) dissipation 
rate for these wavelengths to be smaller than the growth 
rate, which is itself small compared to the shear frequency. 
Therefore, very accurate (e.g. spectral) numerical methods 
are preferable to study this kind of weak instability. 

To check for resolution artefacts, the non stratified in- 
compressible simulations and the compressible simulations 
carried out with NIRVANA below have been computed with 
two different resolutions (the "high" resolution being twice 
the "low" resolution in each direction) . No significant differ- 
ence was found between high resolution and low resolution 
runs. In particular, the same growth rates were obtained, 
with the same localisation properties for the instabilities. 
However, in low resolution runs carried out with NIRVANA 
that were seeded with truncation errors, the unstable scales 
were initially a factor ~ 2 larger, but with the growth rates 
and later nonlinear behaviour being very similar. This can 
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Fig. 7. Vortex streamlines in a sheared flow with x = 3. 
This solution is computed with an FFT method, assuming 
periodic boundary conditions in x and y. This solution and 
the infinite solution (Fig. [1} are extremely similar, except 
around [x — ±5, y — 0) where an X point is observed. 



be understood as a consequence of the scale free property of 
the linear instability that holds once scales are sufficiently 
small but still larger than the dissipation scale, here ex- 
pected to be the grid scale. Noting that conclusions would 
be unaltered if the low resolution runs were adopted, only 
high resolution results will be presented in this section. 

5.2. Numerical solution for an isolated vortex 

To check the existence of the elliptical instability in ac- 
cretion disc vortices, one needs a solution for an isolated 
vortex. In section [51 we have derived an exact solution for a 
vortex in an infinite box. Numerically however, one has to 
work in a finite size domai n, using essentially t he shearing- 
sheet boundary conditions (Ha wlev et al.lll995T h As pointed 
out previously, the infinite solution has a slowly decaying 
tail which can lead to artefacts when one starts with this so- 
lution in a finite-size shearing box. When one tries this pro- 
cedure, the 2D solution rapidly develops instabilities near 
the boundaries, especially along the x axis, leading to large 
amplitude oscillations. To prevent this, we have chosen to 
solve ([8]) using fast Fourier transforms (FFT), which are 
a natural choice for our boundary conditions. A solution 
obtained by this method in a (10 x 4) box is shown in 
Fig. [7] The solution found by this procedure is very sim- 
ilar to the infinite solution, except near (x = ±5, y = 0) 
where an X point is observed. This difference probably ex- 
plains the unsteadiness observed when one uses the infinite 
solution numerically. We find that periodic solutions de- 
rived using the FFT procedure are better suited for numer- 
ical simulations. Although these solutions are not exactly 
steady (mainly because of time varying boundary condi- 
tions), they keep the same shape and the same amplitude 
for several turnover times when one uses high precision nu- 
merical methods, which is enough for our purpose. 

As pointed out previously, this vortex model is not a 
proper infinite elliptical flow, and the approximation used 
in our linear analysis may not be valid. Nevertheless, the 
core of these elliptical vortices is made of elliptical stream- 
lines. Since our linear analysis is essentially a local analysis 
of stability on individual streamlines, one expects our re- 
sults to apply to the core of such vortices. Note however 
that other parametric instabilities may appear outside of 
the vortex core because closed streamlines still exist in this 
region (see Fig. [7]). 

We have computed the evolution of several isolated el- 
liptical vortices in a shearing box with dimensions L x x 
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Fig. 8. Time evolution of y(wf) for several isolated vor- 
tices in shearing box simulations with no stratification. We 
find an exponential growth in agreement with the linear 
prediction for x — 3 and x = 11. We also find an instability 
for x = 5.5 although the elliptical streamlines are stable in 
our linear analysis. 
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Fig. 11. The root mean square vertical velocity obtained 
from the compressible simulations is plotted as a function 
of time for x = 2. The cases shown are v s /c s = 1.3 (dotted 
line), v s /c s = 0.65 (dashed line), and v s /c s — 0.13 (full 
line). Initial values increase with v s /c because the initially 
imposed state variables are further from being in a station- 
ary state. 



L v x L z using an inco mpressible spectral method (see 
iLesur fe Longarettill2005l . for a complete description of the 
numerical method). We have used (N x x N y x N z ) — 
(256 x 512 x 64) points and an eighth order hyperviscosity 
instead of the classical viscosity to prevent the dissipation 
of the weaker unstable modes. The vortices are initialised in 
the center of the box setting a vorticity patch with dimen- 
sion (Aa;, Aj,). We always set X y = 1, L y = 4 and L z = 2, X x 
and L x being adjusted according to the vortex aspect-ratio 
X- For x = 3, we set L x — 10 whereas for x = 5.5 and 
X = 11 we set L x = 20. For each run, we initialise the ve- 
locity field solving §8§ in Fourier space, and we add a small 
(10 -7 ) 3D white noise to the 2D solution. 

5.3. Non stratified case 

To follow the evolution of 3d motions, we have computed 
\J (uf ), where () denotes a volume average procedure. Time 
histories of this quantity are given in Fig. [5] for x = 3, 
X = 5.5 and x = H- As expected from our linear analysis, 
we find an exponential growth for both x = 3 and x = H- 
Since the growth rate for x = 3 is very fast, the instability 
saturates at t ~ 40 and the vortex structure is destroyed, 
as shown in Fig. [9l In the exponential regime, we find a 
growth rate 7 = 0.37 for x = 3 and 7 = 0.016 for x = 11, in 
agreement with theoretical values (73 = 0.44 ; 711 = 0.017). 
We also find an instability for x — 5-5 with 75.5 = 0.055 
which was not expected in our linear analysis (x = 5.5 
elliptical streamlines are stable according to Fig. 2]). 

To check the localisation of the instability, we present 
(x,y) slices of the vertical velocity for x = 3, x = 5.5 and 
X = 1 1 in Fig. |T0j We find that unstable modes are localised 
inside the vortex for x = 3 and x — H, as expected from 
our linear analysis. However, for x — 5.5, the instability 
appears outside of the vortex core, in a region in which one 
find closed non-elliptical streamlines (see Fig. [I}. In this 
particular case, no vertical motions are found in the vortex 
core, which is consistent with our analysis. 

The existence of this outer instability can be understood 
quite simply following the physical description used for the 



elliptical instability fsection l3T4"j) . We know that these insta- 
bilities arise when the epicyclic frequency is an harmonic of 
a closed streamline frequency multiplied by 27r (the stream- 
line frequency being the inverse of the time needed by a 
fluid particle to cover the whole streamline and get back to 
its starting point). Moreover, as one moves away from the 
vortex core, the streamline frequency tends to since the 
flow tends towards a pure shear flow, whereas the epicyclic 
frequency tends to the Keplerian frequency. Therefore, one 
will always be able to find a resonance outside of the vortex 
core, and consequently a parametric instability. 

6. Effects of compressibility 

We have also investigated the effects of compressibil- 
ity. To do this we performed three dimensional hydrody - 
namic simulations using NIRVANA (IZieder fc Yorkd[i~997h . 
NIRVANA has been used frequently in the past to study 
vario us problems involving MHD turbulence in the shearin g 
box (jFromang fc Papaloizoull2006l : iPapaloizou et al.ll2004D . 
Because of its diffusive character, it is best suited to vortices 
with small x that are not too elongated. For this reason we 
limit consideration to x = 2. We consider isothermal shear- 
ing boxes with vertical gravity. In units such that the full 
radial width of the elliptical vorticity patch is unity, the 
box dimensions were (L x — 3.46, L y — 1.73, L z = 2.31). 
The numbers of grid points used were (N x — 196, N y = 
98, N z — 128). A vertical gravitational acceleration — fi 2 z 
was applied over 80% of the vertical computational domain, 
being set to zero in two domains extending from the verti- 
cal boundaries, each being of an extent equal to 10% of the 
whole vertical domain. The boundary conditions of period- 
icity in shearing coordinates could then be applied. We com- 
ment that the use of periodic boundary conditions in the 
vertical direction, with small buffer zones near the bound- 
aries that allow these conditions to be applied consistently, 
has been found to be effective in vertically stra tified shear- 
ing b ox simulations of MHD turbulence (eg. IStone et al.l 
1996). This approach allows for regular treatment of the 
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Fig. 9. 3D snapshots showing the evolution of a \ — 3 vortex at t = 26 (top left), t = 36 (top right), t = 41 (bottom 
left) and t = 47 (bottom right). We have represented in transparent grey an isocontour of vorticity delimiting the vortex 
core, and in colour the volume rendered vertical velocity normalised. We find that the instability is localized inside the 
vortex and is dominated by small scale vertical wave-numbers (0 ~ 0), as expected. Moreover, the vortex structure is 
destroyed when the instability saturates, i.e. for t ~ 40. 



boundary without significant artefacts being produced in 
the mid-plane regions of the flow, as would be expected for 
local instabilities of the type considered here. To test this 
as well as the effects of expanding the computational do- 
main in general, we have performed a test simulation with 
the compuational domain doubled in size in each direction 
for the case with the largest value of v s /c s considered be- 
low. The dimensions of the elliptical vorticity patch were 
not changed. The growth rate and nonlinear development 
of the instability were indeed essentially the same in the 
mid-plane regions, which because of higher inertia carried 
most of the energy, as in the smaller domain simulations 
presented below. The instability in the larger domain case 
readily spreads to the more extensive regions outside the 
original vorticity patch. 

We considered three values of the isothermal sound 
speed c s . Adopting the shearing time as the unit of time, 
these were such that v s /c s were 1.3,0.65, and 0.13. Here 
the velocity v a , being unity in our dimensionless units is 
the maximum velocity difference across the vorticity patch 
due to the background shear. 

The incompressible two dimensional flow field calcu- 
lated in section 15.21 was applied on horizontal planes. As 
above this was calculated by solving the vorticity equa- 
tion with the imposition of periodic boundary conditions. 
To make this possible, one must subtract out the mean 
vorticity in the patch so that the net surface integral 



over the (x,y) plane is zero, which results in some un- 
steadiness as described above. However, in the compress- 
ible case, additional unsteadiness occurs because an asso- 
ciated unbalanced density change is generated. This was 
calculated and included in the set up by noting that for 
a strictly steady state horizontal flow of the type consid- 
ered here, Q(ip) = c 2 hxp+ (l/2)u 2 - (3/2)fi 2 y 2 + $ G , with 
$G = f2 2 z 2 /2 being the gravitational potential, should be 
a function only of the stream function ip. This is readily de- 
termined from the imposed steady flow from the condition 



dQ 1_ , 



(43) 



where u v is taken to be the same function of ip as for the 
infinite vortex patch. The latter is of course an approxima- 
tion as we modified the infinite patch by imposing periodic 
boundary conditions. Nonetheless its use to calculate the 
initial density on horizontal plane resulted in initial states 
that became increasingly more steady as c s increased. 

In each of the simulations the two dimensional vortex 
structure was found to remain for about 40 shearing times 
before undergoing a linear instability which began with a 
very short wavelength in the vertical direction. The root 
mean square (weighted with density) vertical velocity is 
plotted in Fig. 1111 Initial values increase with v s /c because 
the initially imposed flow is further from being stationary. 
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Fig. 10. 2D snapshots of the vertical velocity field in the (x,y) plane for \ — 3 (top), x = 5.5 (middle) and x = 11 
(bottom). The vortex core is delimited by a thick grey line. We find that 3D perturbations are localized inside the vortex 
core for \ — 3 and \ = 11 whereas they are found outside of the vortex core for \ = 5.5. 



Streamlines in the horizontal midplane are shown for 
v s /c s = 1.3 at t = 37, and v s /c s = 0.13 at t = 47 in Fig. [HI 
Although linear instability has begun, the initially imposed 
flow is found to be well preserved for the case with v s /c s = 
0.13 which most closely resembles the incompressible case. 
In this case the maximum growth rate in the linear phase 
of 0.75S' is in good agreement with the linear prediction 
obtained from equation ([2"3"| . 

However, when v s /c s = 1.3 the flow shows significantly 
greater time variability. The vortex is squeezed in the ra- 
dial direction and trailing features begin to form outside 
the initial vortex core which turn into density waves which, 
however, are not well represented in our small computa- 
tional domain. Global vortical structures survive until they 
are eventually destroyed by the instability Snapshots of the 
vertical velocity field in the horizontal midplane are shown 
after the onset of linear instability in Fig. [13] for the three 
simulations. The 3D perturbations are found to be local- 
ized inside the original elliptical boundary of the vorticity 
patch just as in the incompressible case. 

6.1. Stratified case 

To study the vertically stratified case, we have imple- 
mented the Boussinesq equations (|38[) in our spectral code. 
Compared to ([38)) . we have added an eighth order hyper 
diffusivity to the velocity and potential temperature equa- 



tions. We followed the same procedure as in the non strat- 
ified case to set up the vortex structure. We have run the 
simulations with a moderate stratification (N = S), for 
which one expects a departure from the non stratified case 
when x > 4 (Fig. [5]). The time history of (v 2 ) 1 / 2 , for various 
vortex aspect-ratios, is presented in Fig. [T4l 

For x = 3 we find 7 = 0.35 and 7 = 0.01 for x = 5.5. 
No instability is observed in our simulations for % = 11. 
These measured growth rates are always below the ex- 
pected growth rate from our linear analysis (73 = 0.43, 
75.5 = 0.024 and 711 = 0.004) which might be due to exces- 
sive numerical dissipation in the stratified case, despite the 
high resolution and the hyper viscosity prescription used. 
For both 7 = 3 and 7 = 5.5, the instability appears to be 
localised inside the vortex core, as expected. We can con- 
clude from our results that the elliptical instability can be 
detected in stratified simulations, but resolution and nu- 
merical dissipation must be carefully controlled to get ac- 
curate results. 

7. Discussion 

We have described a 3D instability appearing in elliptical 
vortices (Kida vortices) embedded in accretion discs, known 
as the elliptical instability. We have shown analytically, for 
the case of no stratification that this elliptical instability 
is always found in vortex cores, except for a narrow range 
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Fig. 12. Streamlines in the (x, y) plane for x = 2 obtained from the compressible simulations shown after the onset of 
linear instability. These are with v s /c s — 1.3 at t = 37 (left), and v s /c s — 0.13 at t = 47 (right). The initially imposed 
flow is well preserved when v s /c s = 0.13. When v s /c s = 1.3 the flow shows greater time variability, with the vortex being 
sqeezed in the radial direction. Some trailing features are established outside the initial vortex core. 
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Fig. 14. Time evolution of y/ (vf) for several isolated vor- 
tices in stratified shearing box simulations with N = S. We 
find an exponential growth in agreement with the linear 
prediction for \ = 3 and x = 5.5. The expected instability 
for x = 11 is not observed, probably because of a too large 
numerical dissipation. 



of vortex aspect-ratio (4 < x < 5.9). When the vortex is 
weak (elongated vortex), the instability involves small ra- 
dial wavelengths compared to azimuthal and vertical ones, 
the ratio being of the order of x- Moreover, the inclusion 
of a stable stratification tends to weaken the instability for 
large x but it does not suppress it. 

Numerical simulations have essentially confirmed our 
linear analysis. We have found the predicted growth rate in 
the case where the instability was expected inside the core. 
Furthermore, instabilities outside of the vortex core were 
observed when the core was linearly stable, in agreement 
with our physical interpretation. The inclusion of compress- 
ibility in the simulations didn't change the behaviour of the 
instability, and analytical results were recovered. We also 
studied numerically the effect of vertical Boussinesq stratifi- 
cation. Again, analytical results were recovered, except for 



weak vortices (x — 11) for which stability was observed, 
probably because of numerical diffusion. 

These results tend to indicate that elliptical instabili- 
ties, or more generally parametric instabilities, are always 
present in accretion disc vortices. However, they often in- 
volve small radial wavelengths (compared to the disc thick- 
ness) and small growth rates (compared to the orbital time 
scale). Put together, these properties make elliptical insta- 
bilities hard to capture numerically because of numerical 
diffusion at the grid scale. For example, highly accurate 
spectral methods with high-order hyper-diffusivity were re- 
quired to marginally resolve the instability in a x = 11 
vortex. 

In this work, we have not addressed the saturation of 
this instability. This is however rather intentional, as we 
think that the non-linear outcome should be studied to- 
gether with the production mechanism for these vortices. 
In fact, the evolution of these structures will be dictated ul- 
timately by the competition between the production mech- 
anism and the saturated 3D instability. In this context, 
studying the saturation of this instability on its own will 
be of little interest, since it will be modified by the vortex 
production process. Note too that the time-scale needed by 
the saturated instability to destroy the vortex structure is 
unknown. In particular, there is no reason to think that 
this time-scale is related to the growth rate in the linear 
stage. For example, a weakly unstable vortex could be de- 
stroyed in one orbit once saturation is reached. Therefore, 
no firm conclusion can be drawn from this work concern- 
ing the evolution of the vortex structure itself. However, in 
the absence of a sustaining production process, eventually 
one would expect the destruction of the vortex structure 
when the perturbation amplitude reaches the background 
flow amplitude as was indeed observed for our x = 3 simu- 
lation. 

This work has also several limitations. In particular, we 
have assumed an isolated 2D vortex as a starting point. 
However, one might want to study 3D vortices includ- 
ing a complicated vertical structure, with e.g. a mean 
vertical flow in the core (similar to cyclones on Earth 
for example). This kind of structure was suggested by 



12 



G. Lesur and J. C. B. Papaloizou: Stability of Elliptical Vortices 






- 1 .9E-C3 - - 



I 



I 



Fig. 13. 2D snapshots of the vertical velocity field in the (x, y) plane for \ — 2 obtained from the compressible simulations 
shown after the onset of linear instability. These are with v s /c s = 1.3 at t — 37 (top), v s /c s — 0.65 at t = 42 (middle) and 
v s /c s = 0.13 at t = 47. The elliptical boundary of the vorticity patch is overplotted. We find that the 3D perturbations 
are localized inside this boundary as in the incompressible case. 



iBarranco fc Marcus! ((2005) for off-midplane vortices, and 
there is no doubt that stability properties will be modi- 
fied in this case. Moreover, we have mostly restricted our 
study to elliptical streamlines, since these are tractable an- 
alytically. Although some extensions to more general cases 
were shown, we have not provided any formal proof of in- 
stability in the most general case. However, the instability 
mechanism exhibited in the elliptical case appears to be 
quite general, and we think it's unlikely that a steady vor- 
tex could be stable everywhere to 3D perturbations. 

Finally, these results point out the necessity of car- 
rying high-resolution 3D simulations of the mechanisms 
suggested in the literature for vortex formati on. This 
inclu des for example Rossby wave instabilities dLi et ahl 
20011). baroclin i c insta bilities (|Klahr fc Bodenheimerl feOOS: 



Petersen et ail 120071) or gaps due to giant planets 



( de Val-Borro et al.ll2007l ). In any case, the presence of 3D 



parametric instabilities will certainly modify the global out- 
come of these processes, with probably new physical effects. 
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